Approximating real stability radii

The stability radius of an n×n matrix A (or distance to instability) is a well-known measure of robustness of stability of the linear stable dynamical system ẋ = Ax. Such a distance is commonly measured either in the 2-norm or in the Frobenius norm. Even if the matrix A is real, the distance to instability is most often considered with respect to complex valued matrices (in such case the two norms turn out to be equivalent) and restricting the distance to real matrices makes the problem more complicated, and in the case of Frobenius norm to our knowledge unresolved. Here we present a novel approach to approximate real stability radii, particularly well-suited for large sparse matrices. The method consists of a two level iteration, the inner one aiming to compute the εpseudospectral abscissa of a low-rank (1 or 2) dynamical system, and the outer one consisting of an exact Newton iteration. Due to its local convergence property it generally provides upper bounds for the stability radii but in practice usually computes the correct values. The method requires the computation of the rightmost eigenvalue of a sequence of matrices, each of them given by the sum of the original matrix A and a low-rank one. This makes it particularly suitable for large sparse problems, for which several existing methods become inefficient, due to the fact that they require to solve full Hamiltonian eigenvalue problems and/or compute multiple SVDs.